decontam
library(dada2); packageVersion("dada2")## [1] '1.18.0'
library(ggplot2); packageVersion("ggplot2")## [1] '3.3.3'
library(phyloseq); packageVersion("phyloseq")## [1] '1.34.0'
library(data.table); packageVersion("data.table")## [1] '1.13.6'
library(dplyr)
library(plyr)
library(magrittr)
library(RColorBrewer)
library(picante)
library(igraph)
library(cowplot)
library(plotly)
library(lemon)
library(GGally)
library(randomcoloR)
library(knitr)
library(kableExtra)
library(ggpmisc)
library(ggpubr)
library(decontam)
library(stringr)
library(tibble)SVt = readRDS("../seqtab_final.rds")
# SVt_cnn = collapseNoMismatch(SVt, minOverlap = 20,
# identicalOnly = FALSE, vec = TRUE, verbose = TRUE)
# saveRDS(SVt_cnn, "../seqtab_final_cnn.rds")
# SVt = readRDS("../seqtab_final_cnn.rds")
SVt = SVt[!rownames(SVt) %in% "Undetermined", ]
Tax = readRDS("../tax_final.rds")
metadata<-read.csv("../../220319_for_graph.csv", header=TRUE)
# metadata[, 10:47] <- as.numeric(metadata[, 10:47])
metadata$Snowpack_Layer = gsub("Glacial Ice", "GL ICE", metadata$Snowpack_Layer)
metadata$Snowpack_Layer = gsub("Sup. Ice", "SUP ICE", metadata$Snowpack_Layer)
metadata$Snowpack_Layer = gsub("Top", "TOP", metadata$Snowpack_Layer)
metadata$Snowpack_Layer = gsub("Mid", "MID", metadata$Snowpack_Layer)
metadata$Snowpack_Layer = gsub("Basal", "BASAL", metadata$Snowpack_Layer)
order<-rownames(SVt)
metadata$Sample_ID=as.character(metadata$Sample_ID)
metadata_ord<-metadata[match(order, metadata$Sample_ID),]
rownames(metadata_ord) <- metadata_ord$Sample_ID
metadata_ord$Sample_Name <- factor(metadata$Sample_Name, levels=unique(metadata$Sample_Name))
ps <- phyloseq(otu_table(SVt, taxa_are_rows=FALSE),
sample_data(metadata_ord),
tax_table(Tax))
ps.b <- prune_taxa(taxa_sums(ps) > 0, ps)
ps.b = subset_taxa(ps.b, Kingdom %in% c("Archaea", "Bacteria"))
ps.b = subset_samples(ps.b, Fox_Aspect != "N_NE_D")#remove samples with 0 reads
ps.batch = prune_samples(sample_sums(ps.b)>=1, ps.b)
readsumsdf_samples<-data.frame(nreads = sample_sums(ps.batch),
samples = rownames(as.data.frame(sample_sums(ps.batch))),
batch = sample_data(ps.batch)$Batch.no.)
readsumsdf_samples_f = readsumsdf_samples[- grep("Control", readsumsdf_samples$batch),]
my_comparisons=list(c("1","2"),c("1","3"),c("1","4"),c("1","5"),
c("2","3"), c("2","4"), c("2","5"),
c("3","4"), c("3","5"),
c("4","5"))
ggplot(readsumsdf_samples_f, aes(x=batch, y=nreads,
color=batch, fill=batch)) +
geom_violin(alpha=0.3) +
geom_jitter(shape=16, position = "jitter") +
stat_compare_means(comparisons = my_comparisons,
method = "t.test",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
labs(color="Batch Number") +
guides(fill=FALSE) +
xlab("Extraction batch") +
ylab("Number of reads obtained") +
theme(axis.text.x = element_text(vjust = 1,
hjust = 1,
size = 12),
strip.text = element_text(size = 12, margin = margin()))sort(sample_sums(ps.b))## C05-cDNA-29 D12-cDNA-MM C10-cDNA-34 A02-cDNA-MM-H20 C08-cDNA-32
## 0 0 5 7 25
## C07-cDNA-31 B07-cDNA-17 E01-cDNA-MM-H20 B11-cDNA-22 B09-cDNA-19
## 27 43 55 123 187
## D01-cDNA-37 A01-cDNA-MM A11-cDNA-9 D10-cDNA-46 B10-cDNA-21
## 262 290 2320 3844 9170
## A10-cDNA-8 B05-cDNA-15 B08-DNA-4 B11-DNA-7 B08-cDNA-18
## 15595 15781 16781 29350 29808
## C03-cDNA-27 A08-cDNA-6 A06-cDNA-4 A09-cDNA-7 A07-cDNA-5
## 31806 33041 37485 40604 44392
## A11-DNA-9 A05-cDNA-3 A08-DNA-6 B02-cDNA-12 D08-DNA-8
## 50662 52562 53872 57398 57493
## E02-DNA-MM C02-cDNA-26 D07-cDNA-43 D06-DNA-6 D02-DNA-2
## 58537 59570 59805 60653 61230
## A01-DNA-MM B12-DNA-8 B12-cDNA-23 A12-DNA-10 C07-DNA-5
## 61807 62146 62235 63891 64148
## B04-cDNA-14 A07-DNA-5 A02-DNA-MM-H20 A04-cDNA-2 B03-cDNA-13
## 68029 69142 69499 70613 76895
## D09-DNA-9 B02-DNA-12 A05-DNA-3 C04-cDNA-28 C05-DNA-3
## 80028 83134 84732 86983 87568
## C04-DNA-2 D06-cDNA-42 A09-DNA-7 B01-cDNA-11 B04-DNA-14
## 96798 101335 108915 110464 113222
## C03-DNA-B2-1 B03-DNA-13 C08-DNA-6 A12-cDNA-10 B06-cDNA-16
## 114048 118548 122463 123046 123687
## C12-cDNA-36 D02-cDNA-38 D03-DNA-3 C06-DNA-4 E03-DNA-MM-H20
## 123719 137347 138813 139808 140502
## D07-DNA-7 C12-DNA-10 D04-DNA-4 B09-DNA-5 D08-cDNA-44
## 140668 148596 151009 155211 158902
## A06-DNA-4 C01-DNA-9 B10-DNA-6 A10-DNA-8 C09-DNA-7
## 166098 169126 170463 170498 173577
## B01-DNA-11 E01-DNA-4 C06-cDNA-30 C09-cDNA-33 D03-cDNA-39
## 174897 175874 188016 190708 192008
## A03-cDNA-1 A03-DNA-1 B06-DNA-2 C11-cDNA-35 A04-DNA-2
## 193167 196436 198649 198996 200229
## B07-DNA-3 D04-cDNA-40 C10-DNA-8 D11-cDNA-47 D01-DNA-B4-1
## 204577 204668 212618 220186 238199
## D09-cDNA-45 D11-DNA-2 D05-DNA-5 D10-DNA-B5-1 C01-cDNA-25
## 240324 252560 267031 280585 293592
## D05-cDNA-41 D12-DNA-3
## 302761 312011
df <- as.data.frame(sample_data(ps.b)) # Put sample_data into a ggplot-friendly data.frame
df$LibrarySize <- sample_sums(ps.b)
df <- df[order(df$LibrarySize),]
df$Index <- seq(nrow(df))
ggplot(data=df, aes(x=Index, y=LibrarySize, color=Sample_type)) + geom_point()decontam#remove samples with 0 reads
ps.b = prune_samples(sample_sums(ps.b)>=1, ps.b)
sample_data(ps.b)$is.neg <- sample_data(ps.b)$Sample_type == "CONTROL"
contamdf0.5.prev <- isContaminant(ps.b, method="prevalence", threshold = 0.1,
neg="is.neg")
ps.final.noncontam <- prune_taxa(!contamdf0.5.prev$contaminant, ps.b)
ps.final.noncontam.noconts = subset_samples(ps.final.noncontam, Sample_type != "CONTROL")
ps.final.noncontam.noconts.nodoubs = subset_samples(ps.final.noncontam.noconts,
Sample_ID != "D02-DNA-2" & Sample_ID != "B11-DNA-7" & Sample_ID != "C10-cDNA-34" & Sample_ID != "A08-cDNA-6")
ps.final.noncontam.noconts.nodoubs = prune_samples(sample_sums(ps.final.noncontam.noconts.nodoubs)>=1000,
ps.final.noncontam.noconts.nodoubs)
ps.final <- prune_taxa(taxa_sums(ps.final.noncontam.noconts.nodoubs) > 0, ps.final.noncontam.noconts.nodoubs)
ps.final = prune_samples(sample_sums(ps.final)>=1, ps.final)
ps.final.r <- transform_sample_counts(ps.final, function(x) x / sum(x))
ps.final.r.f = filter_taxa(ps.final.r, function(x) sum(x) > .001, TRUE)
ps.final.r## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 35534 taxa and 65 samples ]
## sample_data() Sample Data: [ 65 samples by 47 sample variables ]
## tax_table() Taxonomy Table: [ 35534 taxa by 6 taxonomic ranks ]
ps.final.r.f## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 2886 taxa and 65 samples ]
## sample_data() Sample Data: [ 65 samples by 47 sample variables ]
## tax_table() Taxonomy Table: [ 2886 taxa by 6 taxonomic ranks ]
ps.final## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 35534 taxa and 65 samples ]
## sample_data() Sample Data: [ 65 samples by 47 sample variables ]
## tax_table() Taxonomy Table: [ 35534 taxa by 6 taxonomic ranks ]
sum(rowSums(otu_table(ps.final)))## [1] 6750655
gl_ice_spls = as.character(unique(sample_data(ps.final)[str_detect(sample_data(ps.final)$Sample_Name, "Gl ice"),]$Sample_Name))sample_data(ps.final.r)$Month = factor(sample_data(ps.final.r)$Month,
levels = c("April","June","Early July", "Late July"))
sample_data(ps.final.r)$Fox_Aspect = factor(sample_data(ps.final.r)$Fox_Aspect,
levels = c("NW_AWS","N_NE","SWS1SE"))
sample_data(ps.final.r)$Snowpack_Layer = factor(sample_data(ps.final.r)$Snowpack_Layer,
levels = c("TOP","MID","BASAL","SUP ICE","GL ICE"))
sample_data(ps.final)$Month = factor(sample_data(ps.final)$Month,
levels = c("April","June","Early July", "Late July"))
sample_data(ps.final)$Fox_Aspect = factor(sample_data(ps.final)$Fox_Aspect,
levels = c("NW_AWS","N_NE","SWS1SE"))
sample_data(ps.final)$Snowpack_Layer = factor(sample_data(ps.final)$Snowpack_Layer,
levels = c("TOP","MID","BASAL","SUP ICE","GL ICE"))
sample_data(ps.final.r.f)$Month = factor(sample_data(ps.final.r.f)$Month,
levels = c("April","June","Early July", "Late July"))
sample_data(ps.final.r.f)$Fox_Aspect = factor(sample_data(ps.final.r.f)$Fox_Aspect,
levels = c("NW_AWS","N_NE","SWS1SE"))
sample_data(ps.final.r.f)$Snowpack_Layer = factor(sample_data(ps.final.r.f)$Snowpack_Layer,
levels = c("TOP","MID","BASAL","SUP ICE","GL ICE"))ps.final.r.glom <- tax_glom(subset_samples(ps.final.r,
Snowpack_Layer != "GL ICE" &
Sample_Name != "30 July_SWS1SE_Slush_R" &
Sample_Name != "30 July_SWS1SE_Slush_D" &
!(Sample_Name %in% gl_ice_spls)), taxrank = 'Phylum', NArm = FALSE)
ps.final.r.glom.psdf <- data.table(psmelt(ps.final.r.glom))
ps.final.r.glom.psdf$Phylum <- as.character(ps.final.r.glom.psdf$Phylum)
ps.final.r.glom.psdf[, mean := mean(Abundance, na.rm = TRUE),
by = "Phylum"]
ps.final.r.glom.psdf[(mean <= 0.01), Phylum := "Other"]
ps.final.r.glom.psdf$Phylum[is.na(ps.final.r.glom.psdf$Phylum) & ps.final.r.glom.psdf$Kingdom=="Bacteria"] <- "Unc. Bacteria"
ps.final.r.glom.psdf$Phylum[is.na(ps.final.r.glom.psdf$Phylum) & ps.final.r.glom.psdf$Kingdom=="Archaea"] <- "Unc. Archaea"
#Removing Proteobacteria from previous table
ps.final.r.glom.psdf.noProteo<-ps.final.r.glom.psdf[!grepl("Proteobacteria", ps.final.r.glom.psdf$Phylum),]
#New table, with Proteobacterial classes only
ps.final.r.ProtCl = subset_taxa(ps.final.r, Phylum == "Proteobacteria")
ps.final.r.ProtClglom <- tax_glom(ps.final.r.ProtCl, taxrank = 'Class', NArm = FALSE)
ps.final.r.ProtClglom.psdf <- data.table(psmelt(ps.final.r.ProtClglom))
ps.final.r.ProtClglom.psdf$Class <- as.character(ps.final.r.ProtClglom.psdf$Class)
ps.final.r.ProtClglom.psdf[, mean := mean(Abundance, na.rm = TRUE),
by = "Class"]
ps.final.r.ProtClglom.psdf[(mean <= 0.001), Class := "Other Proteobacteria"]
ps.final.r.ProtClglom.psdf.noP <- subset(ps.final.r.ProtClglom.psdf, select = -Phylum)
names(ps.final.r.glom.psdf.noProteo)[names(ps.final.r.glom.psdf.noProteo) == 'Phylum'] <- 'Taxa'
names(ps.final.r.ProtClglom.psdf.noP)[names(ps.final.r.ProtClglom.psdf.noP) == 'Class'] <- 'Taxa'
# ps.final.r.glom.psdf.noProteo = subset(ps.final.r.glom.psdf.noProteo, select=-c(Class))
ps.final.r.ProteoClAllPhy <- rbind(ps.final.r.glom.psdf.noProteo, ps.final.r.ProtClglom.psdf.noP)
tol18rainbow=c("#771155", "#AA4488", "#CC99BB", "#114477", "#4477AA", "#77AADD", "#117777", "#44AAAA", "#77CCCC", "#777711", "#AAAA44", "#DDDD77", "#774411", "#AA7744", "#DDAA77", "#771122", "#AA4455")
ps.final.r.ProteoClAllPhy$Taxa = factor(ps.final.r.ProteoClAllPhy$Taxa,
levels = c("Alphaproteobacteria", "Deltaproteobacteria",
"Gammaproteobacteria",
"Other Proteobacteria",
"Actinobacteria",
"Acidobacteria","Bacteroidetes",
"Chloroflexi","Cyanobacteria",
"Firmicutes","Gemmatimonadetes","Other"))ggplot(ps.final.r.ProteoClAllPhy,
aes(x = Month, y = Abundance, fill = Taxa)) +
facet_grid(DNA_or_cDNA~Fox_Aspect,
scales = "free",
space = "free") +
geom_bar(stat = "identity",
position = "fill") +
scale_fill_manual(values = tol18rainbow,
na.value = "gray40") +
scale_y_continuous(expand = c(0,0),
labels = scales::percent,
breaks = c(0.25,0.5,0.75,1)) +
theme(axis.title.x = element_blank(),
axis.text.x = element_text(angle = 35, hjust = 1),
strip.text.x = element_text(face = "bold",
size = 12),
strip.background = element_rect(size=1),
legend.text = element_text(size= 12)) +
guides(fill = guide_legend(keywidth = 1, keyheight = 1)) +
ylab("Relative Abundance\n")ggplot(ps.final.r.ProteoClAllPhy,
aes(x = Month, y = Abundance, fill = Taxa)) +
facet_grid(DNA_or_cDNA~Snowpack_Layer,
scales = "free",
space = "free") +
geom_bar(stat = "identity",
position = "fill") +
scale_fill_manual(values = tol18rainbow,
na.value = "gray40") +
scale_y_continuous(expand = c(0,0),
labels = scales::percent,
breaks = c(0.25,0.5,0.75,1)) +
theme(axis.title.x = element_blank(),
axis.text.x = element_text(angle = 35, hjust = 1),
strip.text.x = element_text(face = "bold",
size = 11),
strip.background = element_rect(size=1),
legend.text = element_text(size= 12)) +
guides(fill = guide_legend(keywidth = 1, keyheight = 1)) +
ylab("Relative Abundance\n")tol17rainbow=c("#771155", "#AA4488", "#CC99BB", "#114477", "#4477AA", "#77AADD", "#117777", "#44AAAA", "#77CCCC", "#777711", "#AAAA44", "#DDDD77", "#774411", "#AA7744", "#DDAA77", "#AA4455")
ggplot(ps.final.r.ProteoClAllPhy,
aes(x = Sample_Name, y = Abundance, fill = Taxa)) +
facet_grid(DNA_or_cDNA~Month,
scales = "free",
space = "free") +
geom_bar(stat = "identity",
position = "fill",
color="black") +
scale_fill_manual(values = tol17rainbow,
na.value = "gray40") +
scale_y_continuous(expand = c(0,0),
labels = scales::percent,
breaks = c(0.25,0.5,0.75,1)) +
theme(axis.title.x = element_blank(),
axis.text.x = element_text(angle = 35, hjust = 1),
strip.text.x = element_text(face = "bold",
size = 16),
strip.background = element_rect(size=1),
legend.text = element_text(size= 14)) +
guides(fill = guide_legend(keywidth = 1, keyheight = 1,
nrow = )) +
ylab("Relative Abundance\n")Top20SVs_names = names(sort(taxa_sums(ps.final.r), TRUE)[1:20])
Top20SVs_names_df=as.data.frame(as(tax_table(prune_taxa(Top20SVs_names, ps.final.r)), "matrix")[,2:6])
rownames(Top20SVs_names_df) = c()
Top20SVs_names_df## Phylum Class Order Family
## 1 Proteobacteria Alphaproteobacteria Rhizobiales Beijerinckiaceae
## 2 Proteobacteria Gammaproteobacteria Betaproteobacteriales Burkholderiaceae
## 3 Proteobacteria Alphaproteobacteria Rhizobiales Labraceae
## 4 Proteobacteria Alphaproteobacteria Sphingomonadales Sphingomonadaceae
## 5 Proteobacteria Gammaproteobacteria Betaproteobacteriales Burkholderiaceae
## 6 Proteobacteria Alphaproteobacteria Rhizobiales Xanthobacteraceae
## 7 Proteobacteria Alphaproteobacteria Sphingomonadales Sphingomonadaceae
## 8 Proteobacteria Alphaproteobacteria Rhizobiales Beijerinckiaceae
## 9 Proteobacteria Alphaproteobacteria Sphingomonadales Sphingomonadaceae
## 10 Proteobacteria Alphaproteobacteria Acetobacterales Acetobacteraceae
## 11 Actinobacteria Actinobacteria Micrococcales Microbacteriaceae
## 12 Proteobacteria Gammaproteobacteria Betaproteobacteriales Burkholderiaceae
## 13 Proteobacteria Alphaproteobacteria Sphingomonadales Sphingomonadaceae
## 14 Proteobacteria Gammaproteobacteria Betaproteobacteriales Burkholderiaceae
## 15 Actinobacteria Actinobacteria Micrococcales Micrococcaceae
## 16 Actinobacteria Actinobacteria Corynebacteriales Nocardiaceae
## 17 Actinobacteria Actinobacteria Frankiales <NA>
## 18 Actinobacteria Actinobacteria Micrococcales Cellulomonadaceae
## 19 Actinobacteria Actinobacteria Frankiales <NA>
## 20 Actinobacteria Actinobacteria Frankiales <NA>
## Genus
## 1 Bosea
## 2 Cupriavidus
## 3 Labrys
## 4 Sphingomonas
## 5 Burkholderia-Caballeronia-Paraburkholderia
## 6 Bradyrhizobium
## 7 Sphingomonas
## 8 Methylobacterium
## 9 Sphingomonas
## 10 <NA>
## 11 Salinibacterium
## 12 Massilia
## 13 Sphingomonas
## 14 Delftia
## 15 Arthrobacter
## 16 <NA>
## 17 <NA>
## 18 <NA>
## 19 <NA>
## 20 <NA>
rownames(as.data.frame(as(tax_table(prune_taxa(Top20SVs_names, ps.final.r)), "matrix")[,2:6]))## [1] "TGGGGAATATTGGACAATGGGCGCAAGCCTGATCCAGCCATGCCGCGTGAGTGATGAAGGCCTTAGGGTTGTAAAGCTCTTTTGTCCGGGAAGATAATGACTGTACCGGAAGAATAAGCCCCGGCTAACTTCGTGCCAGCAGCCGCGGTAATACGAAGGGGGCTAGCGTTGCTCGGAATCACTGGGCGTAAAGGGCGCGTAGGCGGACTCTTAAGTCGGGGGTGAAAGCCCAGGGCTCAACCCTGGAATTGCCTTCGATACTGGGAGTCTTGAGTTCGGAAGAGGTTGGTGGAACTGCGAGTGTAGAGGTGAAATTCGTAGATATTCGCAAGAACACCGGTGGCGAAGGCGGCCAACTGGTCCGAAACTGACGCTGAGGCGCGAAAGCGTGGGGAGCAAACA"
## [2] "TGGGGAATTTTGGACAATGGGGGCAACCCTGATCCAGCAATGCCGCGTGTGTGAAGAAGGCCTTCGGGTTGTAAAGCACTTTTGTCCGGAAAGAAATCGCGCTGGTTAATACCTGGCGTGGATGACGGTACCGGAAGAATAAGCACCGGCTAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGTGCGAGCGTTAATCGGAATTACTGGGCGTAAAGCGTGCGCAGGCGGTTTTGTAAGACAGGCGTGAAATCCCCGGGCTTAACCTGGGAATTGCGCTTGTGACTGCAAGGCTAGAGTGCGTCAGAGGGGGGTAGAATTCCACGTGTAGCAGTGAAATGCGTAGAGATGTGGAGGAATACCGATGGCGAAGGCAGCCCCCTGGGACGTGACTGACGCTCATGCACGAAAGCGTGGGGAGCAAACA"
## [3] "TGGGGAATATTGGACAATGGGCGCAAGCCTGATCCAGCCATGCCGCGTGAGTGATGACGGCCTTAGGGTTGTAAAGCTCTTTTAACAGGGACGATAATGACGGTACCTGTAGAATAAGCCCCGGCAAACTTCGTGCCAGCAGCCGCGGTAATACGAAGGGGGCTAGCGTTGTTCGGAATTACTGGGCGTAAAGCGCACGTAGGCGGATTGTTAAGTCGGGGGTGAAATCCTGAGGCTCAACCTCAGAACTGCCTTCGATACTGGCGATCTTGAGTTCGGAAGAGGTTGGTGGAACAGCTAGTGTAGAGGTGAAATTCGTAGATATTAGCTAGAACACCAGTGGCGAAGGCGGCCAACTGGTCCGATACTGACGCTGAGGTGCGAAAGCGTGGGGAGCAAACA"
## [4] "TGGGGAATATTGGACAATGGGCGAAAGCCTGATCCAGCAATGCCGCGTGAGTGATGAAGGCCTTAGGGTTGTAAAGCTCTTTTACCCGGGATGATAATGACAGTACCGGGAGAATAAGCTCCGGCTAACTCCGTGCCAGCAGCCGCGGTAATACGGAGGGAGCTAGCGTTATTCGGAATTACTGGGCGTAAAGCGCACGTAGGCGGCTTTGTAAGTAAGAGGTGAAAGCCTGGTGCTCAACACCAGAACTGCCTTTTAGACTGCATCGCTTGAATCCAGGAGAGGTGAGTGGAATTCCGAGTGTAGAGGTGAAATTCGTAGATATTCGGAAGAACACCAGTGGCGAAGGCGGCTCACTGGACTGGTATTGACGCTGAGGTGCGAAAGCGTGGGGAGCAAACA"
## [5] "TGGGGAATTTTGGACAATGGGGGAAACCCTGATCCAGCAATGCCGCGTGTGTGAAGAAGGCCTTCGGGTTGTAAAGCACTTTTGTCCGGAAAGAAAACCTTTGCACTAATACTGTGAGGGGATGACGGTACCGGAAGAATAAGCACCGGCTAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGTGCGAGCGTTAATCGGAATTACTGGGCGTAAAGCGTGCGCAGGCGGTTCGTTAAGACAGATGTGAAATCCCCGGGCTTAACCTGGGAACTGCATTTGTGACTGGCGAGCTAGAGTATGGCAGAGGGGGGTAGAATTCCACGTGTAGCAGTGAAATGCGTAGAGATGTGGAGGAATACCGATGGCGAAGGCAGCCCCCTGGGCCAATACTGACGCTCATGCACGAAAGCGTGGGGAGCAAACA"
## [6] "TGGGGAATATTGGACAATGGGCGCAAGCCTGATCCAGCCATGCCGCGTGAGTGATGAAGGCCCTAGGGTTGTAAAGCTCTTTTGTGCGGGAAGATAATGACGGTACCGCAAGAATAAGCCCCGGCTAACTTCGTGCCAGCAGCCGCGGTAATACGAAGGGGGCTAGCGTTGCTCGGAATCACTGGGCGTAAAGGGTGCGTAGGCGGGTCTTTAAGTCAGGGGTGAAATCCTGGAGCTCAACTCCAGAACTGCCTTTGATACTGAAGATCTTGAGTTCGGGAGAGGTGAGTGGAACTGCGAGTGTAGAGGTGAAATTCGTAGATATTCGCAAGAACACCAGTGGCGAAGGCGGCTCACTGGCCCGATACTGACGCTGAGGCACGAAAGCGTGGGGAGCAAACA"
## [7] "TGGGGAATATTGGACAATGGGCGAAAGCCTGATCCAGCAATGCCGCGTGAGTGATGAAGGCCTTAGGGTTGTAAAGCTCTTTTACCCGGGATGATAATGACAGTACCGGGAGAATAAGCTCCGGCTAACTCCGTGCCAGCAGCCGCGGTAATACGGAGGGAGCTAGCGTTATTCGGAATTACTGGGCGTAAAGCGCACGTAGGCGGCTTTGTAAGTAAGAGGTGAAAGCCCAGAGCTCAACTCTGGAATTGCCTTTTAGACTGCATCGCTTGAATCATGGAGAGGTCAGTGGAATTCCGAGTGTAGAGGTGAAATTCGTAGATATTCGGAAGAACACCAGTGGCGAAGGCGGCTGACTGGACATGTATTGACGCTGAGGTGCGAAAGCGTGGGGAGCAAACA"
## [8] "TGGGGAATATTGGACAATGGGCGCAAGCCTGATCCAGCCATGCCGCGTGAGTGATGAAGGCCTTAGGGTTGTAAAGCTCTTTTATCCGGGACGATAATGACGGTACCGGAGGAATAAGCCCCGGCTAACTTCGTGCCAGCAGCCGCGGTAATACGAAGGGGGCTAGCGTTGCTCGGAATCACTGGGCGTAAAGGGCGCGTAGGCGGCGTTTTAAGTCGGGGGTGAAAGCCTGTGGCTCAACCACAGAATGGCCTTCGATACTGGGACGCTTGAGTATGGTAGAGGTTGGTGGAACTGCGAGTGTAGAGGTGAAATTCGTAGATATTCGCAAGAACACCGGTGGCGAAGGCGGCCAACTGGACCATCACTGACGCTGAGGCGCGAAAGCGTGGGGAGCAAACA"
## [9] "TGGGGAATATTGGACAATGGGCGAAAGCCTGATCCAGCAATGCCGCGTGAGTGATGAAGGCCCTAGGGTTGTAAAGCTCTTTTACCCGGGAAGATAATGACTGTACCGGGAGAATAAGCCCCGGCTAACTCCGTGCCAGCAGCCGCGGTAATACGGAGGGGGCTAGCGTTGTTCGGAATTACTGGGCGTAAAGCGCACGTAGGCGGCTTTGTAAGTCAGAGGTGAAAGCCTGGAGCTCAACTCCAGAACTGCCTTTGAGACTGCATCGCTTGAATCCAGGAGAGGTCAGTGGAATTCCGAGTGTAGAGGTGAAATTCGTAGATATTCGGAAGAACACCAGTGGCGAAGGCGGCTGACTGGACTGGTATTGACGCTGAGGTGCGAAAGCGTGGGGAGCAAACA"
## [10] "TGGGGAATATTGGACAATGGGCGAAAGCCTGATCCAGCAATGCCGCGTGTGTGAAGAAGGTCTTCGGATTGTAAAGCACTTTCGGCGGGGACGATGATGACGGTACCCGCAGAAGAAGCCCCGGCTAACTTCGTGCCAGCAGCCGCGGTAATACGAAGGGGGCTAGCGTTGCTCGGAATGACTGGGCGTAAAGGGCGCGTAGGCGGTTGCATAAGTTAGATGTGAAATTCCCGGGCTTAACCTGGGGACTGCATTTGATACTATGTGGCTTGAGTATGGAAGAGGGTCGTGGAATTCCCAGTGTAGAGGTGAAATTCGTAGATATTGGGAAGAACACCGGTGGCGAAGGCGGCGACCTGGTCCATAACTGACGCTGAGGCGCGAAAGCGTGGGGAGCAAACA"
## [11] "TGGGGAATATTGCACAATGGGCGCAAGCCTGATGCAGCAACGCCGCGTGAGGGACGACGGCCTTCGGGTTGTAAACCTCTTTTAGTAGGGAAGAAGCGAAAGTGACGGTACCTGCAGAAAAAGCACCGGCTAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGTGCAAGCGTTATCCGGAATTATTGGGCGTAAAGAGCTCGTAGGCGGTTTGTCGCGTCTGCTGTGAAAACTGGGGGCTCAACCCCCAGCCTGCAGTGGGTACGGGCAGACTAGAGTGCGGTAGGGGAGATTGGAATTCCTGGTGTAGCGGTGGAATGCGCAGATATCAGGAGGAACACCAATGGCGAAGGCAGATCTCTGGGCCGTAACTGACGCTGAGGAGCGAAAGCATGGGGAGCGAACA"
## [12] "TGGGGAATTTTGGACAATGGGCGAAAGCCTGATCCAGCAATGCCGCGTGAGTGAAGAAGGCCTTCGGGTTGTAAAGCTCTTTTGTCAGGGAAGAAACGGTGTGGGCTAATATCCTGCACTAATGACGGTACCTGAAGAATAAGCACCGGCTAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGTGCGCAGGCGGTTTTGTAAGTCTGTCGTGAAATCCCCGGGCTCAACCTGGGAATTGCGATGGAGACTGCAAGGCTAGAATCTGGCAGAGGGGGGTAGAATTCCACGTGTAGCAGTGAAATGCGTAGAGATGTGGAGGAACACCGATGGCGAAGGCAGCCCCCTGGGTCAAGATTGACGCTCATGCACGAAAGCGTGGGGAGCAAACA"
## [13] "TGGGGAATATTGGACAATGGGCGAAAGCCTGATCCAGCAATGCCGCGTGAGTGATGAAGGCCTTAGGGTTGTAAAGCTCTTTTACCCGGGATGATAATGACAGTACCGGGAGAATAAGCTCCGGCTAACTCCGTGCCAGCAGCCGCGGTAATACGGAGGGAGCTAGCGTTATTCGGAATTACTGGGCGTAAAGCGCACGTAGGCGGCTTTGTAAGTTAGAGGTGAAAGCCTGGAGCTCAACTCCAGAATTGCCTTTAAGACTGCATCGCTTGAATCCAGGAGAGGTGAGTGGAATTCCGAGTGTAGAGGTGAAATTCGTAGATATTCGGAAGAACACCAGTGGCGAAGGCGGCTCACTGGACTGGTATTGACGCTGAGGTGCGAAAGCGTGGGGAGCAAACA"
## [14] "TGGGGAATTTTGGACAATGGGCGAAAGCCTGATCCAGCAATGCCGCGTGCAGGATGAAGGCCTTCGGGTTGTAAACTGCTTTTGTACGGAACGAAAAAGCTCCTTCTAATACAGGGGGCCCATGACGGTACCGTAAGAATAAGCACCGGCTAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGTGCGAGCGTTAATCGGAATTACTGGGCGTAAAGCGTGCGCAGGCGGTTATGTAAGACAGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATTTGTGACTGCATGGCTAGAGTACGGTAGAGGGGGATGGAATTCCGCGTGTAGCAGTGAAATGCGTAGATATGCGGAGGAACACCGATGGCGAAGGCAATCCCCTGGACCTGTACTGACGCTCATGCACGAAAGCGTGGGGAGCAAACA"
## [15] "TGGGGAATATTGCACAATGGGCGAAAGCCTGATGCAGCGACGCCGCGTGAGGGATGACGGCCTTCGGGTTGTAAACCTCTTTCAGTAGGGAACAAGGCCACACGTGGTGTGGTTGAGGGTACTTGCAGAAGAAGCGCCGGCTAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGCGCAAGCGTTATCCGGAATTATTGGGCGTAAAGAGCTCGTAGGCGGTTTGTCGCGTCTGCCGTGAAAGTCCGGGGCTCAACTCCGGATCTGCGGTGGGTACGGGCAGACTAGAGTGATGTAGGGGAGACTGGAATTCCTGGTGTAGCGGTGAAATGCGCAGATATCAGGAGGAACACCGATGGCGAAGGCAGGTCTCTGGGCATTAACTGACGCTGAGGAGCGAAAGCATGGGGAGCGAACA"
## [16] "TGGGGAATATTGCACAATGGGCGAAAGCCTGATGCAGCGACGCCGCGTGAGGGATGACGGCCTTCGGGTTGTAAACCTCTTTCAGCAGGGACGAAGCGAGAGTGACGGTACCTGCAGAAGAAGCACCGGCCAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGTGCGAGCGTTGTCCGGAATTACTGGGCGTAAAGAGCTCGTAGGCGGTTTGTCACGTCGGCTGTGAAAACCCATCGCTCAACGGTGGGCTTGCAGTCGATACGGGCTGACTTGAGTACTGCAGGGGAGACTGGAATTCCTGGTGTAGCGGTGAAATGCGCAGATATCAGGAGGAACACCGGTGGCGAAGGCGGGTCTCTGGGCAGTAACTGACGCTGAGGAGCGAAAGCGTGGGGAGCGAACA"
## [17] "TGGGGAATATTGCGCAATGGGCGAAAGCCTGACGCAGCGACGCCGCGTGGGGGATGAAGGCCTTCGGGTTGTAAACCTCTTTCAGCTCCGACGAATTCAGACGGTAGGAGCAGAAGAAGCACCGGCCAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGTGCAAGCGTTGTCCGGAATTATTGGGCGTAAAGAGCTCGTAGGCGGTTTGTCGCGTCGGCTGTGAAATCCCGTGGCTCAACTGCGGGTCTGCAGTCGATACGGGCAGACTAGAGTGCGGTAGGGGAGACTGGAATTCCTGGTGTAGCGGTGAAATGCGCAGATATCAGGAGGAACACCGGTGGCGAAGGCGGGTCTCTGGGCCGTAACTGACGCTGAGGAGCGAAAGCGTGGGGAGCGAACA"
## [18] "TGGGGAATATTGCACAATGGGCGAAAGCCTGATGCAGCGACGCCGCGTGAGGGATGGAGGCCTTCGGGTTGTAAACCTCTTTCAGCAGGGAAGAAGTGGTCGGGTTCTCTCGGTCGTGACGGTACCTGCAGAAGAAGCGCCGGCTAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGCGCAAGCGTTGTCCGGAATTATTGGGCGTAAAGAGCTCGTAGGCGGTTTGTCGCGTCTGCTGTGAAAACTCGAGGCTCAACCTCGGGCTTGCAGTGGGTACGGGCAGACTAGAGTGCGGTAGGGGAGACTGGAATTCCTGGTGTAGCGGTGGAATGCGCAGATATCAGGAGGAACACCGATGGCGAAGGCAGGTCTCTGGGCCGCAACTGACGCTGAGGAGCGAAAGCATGGGGAGCGAACA"
## [19] "TGGGGAATATTGCGCAATGGGCGAAAGCCTGACGCAGCGACGCCGCGTGGGGGATGAAGGCCTTCGGGTTGTAAACCTCTTTCAGCTCCGACGAATTCAGACGGTAGGAGCAGAAGAAGCACCGGCCAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGTGCAAGCGTTGTCCGGAATTATTGGGCGTAAAGAGCTCGTAGGCGGTTTGTCGCGTCGGCTGTGAAATCCCGTGGCTCAACTGCGGGTCTGCAGTCGATACGGGCAGACTAGAGTGCGGTAGGGGAGACTGGAATTCCTGGTGTAGCGGTGAAATGCGCAGATATCAGGAGGAACACCGGTGGCGAAGGCGGGTCTCTGGGCCGTAACTGACGCTGAGGAGCGAAAGCGTGGGGAGCAAACA"
## [20] "TGGGGAATATTGCGCAATGGGCGAAAGCCTGACGCAGCGACGCCGCGTGAGGGATGACGGCCTTCGGGTTGTAAACCTCTTTCAGCAGGGACGAACACAGACGGTACCTGCAGAAGAAGCACCGGCCAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGTGCAAGCGTTGTCCGGAATTATTGGGCGTAAAGAGCTCGTAGGCGGTCTGTCGCGTCGGCTGTGAAAACTCGGGGCTCAACCCCGAGCCTGCAGTCGATACGGGCAGACTAGAGTGCTGTAGGGGAGACTGGAATTCCTGGTGTAGCGGTGAAATGCGCAGATATCAGGAGGAACACCGGTGGCGAAGGCGGGTCTCTGGGCAGTAACTGACGCTGAGGAGCGAAAGCGTGGGGAGCGAACA"
Oh good. But this need improving. Getting the unclassified component of the families for the top main taxa in it, but without a guarantee for them popping into it. Also, basing it on the >0.1% dataset as per Arwyn’s advice.
ps.final.r.glom.genus.psdf_2 <- data.table(speedyseq::psmelt(ps.final.r.glom.genus))
ps.final.r.glom.genus.psdf_2$Genus <- as.character(ps.final.r.glom.genus.psdf_2$Genus)
ps.final.r.glom.genus.psdf_2[, mean := mean(Abundance, na.rm = TRUE),
by = "Genus"]
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Sphingomonadaceae"] <- "Unc. Sphingomonadaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Burkholderiaceae"] <- "Unc. Burkholderiaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Labraceae"] <- "Unc. Labraceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Acetobacteraceae"] <- "Unc. Acetobacteraceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Microbacteriaceae"] <- "Unc. Microbacteriaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Deinococcaceae"] <- "Unc. Deinococcaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Gemmatimonaceae"] <- "Unc. Gemmatimonaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Intrasporangiaceae"] <- "Unc. Intrasporangiaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Nocardiaceae"] <- "Unc. Nocardiaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Ktedonobacteraceae"] <- "Unc. Ktedonobacteraceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Cellulomonadaceae"] <- "Unc. Cellulomonadaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & ps.final.r.glom.genus.psdf_2$Family=="Beijerinckiaceae"] <- "Unc. Beijerinckiaceae"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & is.na(ps.final.r.glom.genus.psdf_2$Family)] <- "Unclassified"
ps.final.r.glom.genus.psdf_2$Genus[is.na(ps.final.r.glom.genus.psdf_2$Genus) & is.na(ps.final.r.glom.genus.psdf_2$Family) & ps.final.r.glom.genus.psdf_2$Order=="Frankiales"] <- "Unc. Frankiales"
ps.final.r.glom.genus.psdf_2[(mean <= 0.01), Genus := "Other (<1%)"]
ps.final.r.glom.genus.psdf_2$Genus = factor(ps.final.r.glom.genus.psdf_2$Genus,
levels = c("Acidiphilium","Arthrobacter",
"Bacillus",
"Bosea",
"Burkholderia-Caballeronia-Paraburkholderia",
"Cryobacterium","Cupriavidus",
"Delftia","Hymenobacter",
"Labrys","Leptothrix",
"Massilia","Methylobacterium",
"Noviherbaspirillum","Oryzihumus",
"Polymorphobacter",
"Salinibacterium",
"Sphingomonas","Other (<1%)"))tol17rainbow=c("#771155", "#AA4488", "#CC99BB", "#114477",
"#4477AA", "#77AADD", "#117777", "#44AAAA",
"#77CCCC", "#777711", "#AAAA44", "#DDDD77",
"#774411", "#AA7744", "#DDAA77", "#771122",
"#AA4455", "#DD7788", "azure3")
ggplot(ps.final.r.glom.genus.psdf_2,
aes(x = Month, y = Abundance, fill = Genus)) +
facet_grid(DNA_or_cDNA~Snowpack_Layer,
scales = "free",
space = "free") +
geom_bar(stat = "identity",
position = "fill",
color="black") +
scale_fill_manual(values = tol17rainbow,
na.value = "gray40") +
scale_y_continuous(expand = c(0,0),
labels = scales::percent,
breaks = c(0.25,0.5,0.75,1)) +
theme(axis.title.x = element_blank(),
axis.text.x = element_text(angle = 35, hjust = 1),
strip.text.x = element_text(face = "bold",
size = 16),
strip.background = element_rect(size=1),
legend.text = element_text(size= 14)) +
guides(fill = guide_legend(keywidth = 1, keyheight = 1,
nrow = )) +
ylab("Relative Abundance\n")ggplot(ps.final.r.glom.genus.psdf_2,
aes(x = Month, y = Abundance, fill = Genus)) +
facet_grid(DNA_or_cDNA~Fox_Aspect,
scales = "free",
space = "free") +
geom_bar(stat = "identity",
position = "fill",
color="black") +
scale_fill_manual(values = tol17rainbow,
na.value = "gray40") +
scale_y_continuous(expand = c(0,0),
labels = scales::percent,
breaks = c(0.25,0.5,0.75,1)) +
theme(axis.title.x = element_blank(),
axis.text.x = element_text(angle = 35, hjust = 1),
strip.text.x = element_text(face = "bold",
size = 16),
strip.background = element_rect(size=1),
legend.text = element_text(size= 14)) +
guides(fill = guide_legend(keywidth = 1, keyheight = 1,
nrow = )) +
ylab("Relative Abundance\n")Lest us not forget that size of Other (<1%) corresponds to the different samples included in that category.
I don’t like this, but plotting all samples now, without grouping.
# ps.final.r.glom.genus.psdf_2$Sample_Name <- as.character(ps.final.r.glom.genus.psdf_2$Sample_Name)
# ps.final.r.glom.genus.psdf_2$Sample_Name <- factor(ps.final.r.glom.genus.psdf_2$Sample_Name,
# levels=unique(ps.final.r.glom.genus.psdf_2$Sample_Name))
ggplot(ps.final.r.glom.genus.psdf_2,
aes(x = Sample_Name, y = Abundance, fill = Genus)) +
facet_grid(DNA_or_cDNA~Month,
scales = "free",
space = "free") +
geom_bar(stat = "identity",
position = "fill",
color="black") +
scale_fill_manual(values = tol17rainbow,
na.value = "gray40") +
scale_y_continuous(expand = c(0,0),
labels = scales::percent,
breaks = c(0.25,0.5,0.75,1)) +
theme(axis.title.x = element_blank(),
axis.text.x = element_text(angle = 35, hjust = 1),
strip.text.x = element_text(face = "bold",
size = 16),
strip.background = element_rect(size=1),
legend.text = element_text(size= 14)) +
guides(fill = guide_legend(keywidth = 1, keyheight = 1,
nrow = )) +
ylab("Relative Abundance\n")nmds_shapes = c(15,16,17,18,6)
nmdscolors_by_month = c("April" = "#A6048B",
"June" = "#0704A6",
"Early July" = "#727A0B",
"Late July" = "#6C0504")
nmdscolors_by_snpl = c("Top" = "#151AC6",
"Mid" = "#C6156D",
"Basal" = "#890909",
"Sup. Ice" = "#091D89",
"Glacial Ice" = "#580303")
ps.final.r.dna = subset_samples(ps.final.r, DNA_or_cDNA != "cDNA")
ps.final.r.cdna = subset_samples(ps.final.r, DNA_or_cDNA != "DNA")
ps.final.r.dna.nmds <- ordinate(ps.final.r.dna,
"NMDS",
distance="jsd",
maxit = 1e3)## Run 0 stress 0.2107623
## Run 1 stress 0.2270207
## Run 2 stress 0.2259277
## Run 3 stress 0.2158645
## Run 4 stress 0.2249122
## Run 5 stress 0.2276945
## Run 6 stress 0.2254988
## Run 7 stress 0.2087417
## ... New best solution
## ... Procrustes: rmse 0.1326501 max resid 0.2782618
## Run 8 stress 0.2034057
## ... New best solution
## ... Procrustes: rmse 0.0775235 max resid 0.4113966
## Run 9 stress 0.2030457
## ... New best solution
## ... Procrustes: rmse 0.1023723 max resid 0.4191121
## Run 10 stress 0.2178268
## Run 11 stress 0.2064812
## Run 12 stress 0.2106171
## Run 13 stress 0.2087421
## Run 14 stress 0.2030472
## ... Procrustes: rmse 0.0006936713 max resid 0.00327211
## ... Similar to previous best
## Run 15 stress 0.2264449
## Run 16 stress 0.2019616
## ... New best solution
## ... Procrustes: rmse 0.04318994 max resid 0.188176
## Run 17 stress 0.208532
## Run 18 stress 0.2236278
## Run 19 stress 0.2161764
## Run 20 stress 0.2188811
## *** No convergence -- monoMDS stopping criteria:
## 1: no. of iterations >= maxit
## 18: stress ratio > sratmax
## 1: scale factor of the gradient < sfgrmin
ps.final.r.cdna.nmds <- ordinate(ps.final.r.cdna,
"NMDS",
distance="jsd",
maxit = 1e3)## Run 0 stress 0.1654426
## Run 1 stress 0.1595584
## ... New best solution
## ... Procrustes: rmse 0.1139479 max resid 0.3178704
## Run 2 stress 0.1710025
## Run 3 stress 0.1718853
## Run 4 stress 0.1648546
## Run 5 stress 0.1678885
## Run 6 stress 0.1649862
## Run 7 stress 0.1686493
## Run 8 stress 0.1656653
## Run 9 stress 0.166164
## Run 10 stress 0.1599933
## ... Procrustes: rmse 0.0790213 max resid 0.2159161
## Run 11 stress 0.162065
## Run 12 stress 0.1724002
## Run 13 stress 0.176513
## Run 14 stress 0.1667006
## Run 15 stress 0.1616702
## Run 16 stress 0.1697515
## Run 17 stress 0.1661504
## Run 18 stress 0.1651428
## Run 19 stress 0.1605854
## Run 20 stress 0.1687601
## *** No convergence -- monoMDS stopping criteria:
## 19: stress ratio > sratmax
## 1: scale factor of the gradient < sfgrmin
ps.final.r.dna.nmds.plot <- plot_ordination(ps.final.r.dna, ps.final.r.dna.nmds) +
geom_point(aes(colour=Month, shape=Snowpack_Layer),
size=5, alpha = 0.5, stroke=1.5) +
scale_shape_manual(values=nmds_shapes) +
scale_colour_manual(values=nmdscolors_by_month) +
scale_fill_manual(values=nmdscolors_by_month) +
stat_ellipse(aes(fill=Month),
geom = "polygon",
type="t",
alpha=0.15,
linetype = 2,
show.legend = FALSE) +
ggtitle("DNA nMDS") +
guides(colour = guide_legend(title="Month"))
ps.final.r.cdna.nmds.plot <- plot_ordination(ps.final.r.cdna, ps.final.r.cdna.nmds) +
geom_point(aes(colour=Month, shape=Snowpack_Layer),
size=5, alpha = 0.7, stroke=1.5) +
scale_shape_manual(values=nmds_shapes) +
scale_colour_manual(values=nmdscolors_by_month) +
scale_fill_manual(values=nmdscolors_by_month) +
stat_ellipse(aes(fill=Month),
geom = "polygon",
type="t",
alpha=0.2,
linetype = 2,
show.legend = FALSE)+
ggtitle("cDNA nMDS") +
guides(colour = guide_legend(title="Month"))plot_grid(ps.final.r.dna.nmds.plot + theme_bw(),
ps.final.r.cdna.nmds.plot + theme_bw())month_comparisons = list( c("April", "June"), c("April", "Early July"),
c("April", "Late July"),
c("June", "Early July"), c("June", "Late July"),
c("Early July", "Late July"))
site_comparisons = list( c("NW_AWS", "N_NE"),
c("NW_AWS", "SWS1SE"),
c("N_NE", "SWS1SE"))
slayer_comparisons = list( c("TOP", "MID"), c("TOP", "BASAL"), c("TOP", "SUP ICE"),
c("TOP", "GL ICE"),
c("MID", "BASAL"), c("MID", "SUP ICE"),
c("MID", "GL ICE"),
c("BASAL", "SUP ICE"), c("BASAL", "GL ICE"),
c("SUP ICE", "GL ICE"))plot_richness(ps.final,
x="Sample_ID", color="Sample_ID",
measures=c("Observed","Shannon","Simpson")) +
geom_jitter(shape=16, position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
labs(color="Sample_ID") +
guides(color=guide_legend(ncol=2)) +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))mini_metadata = metadata %>%
select(c(Sample_ID, Sample_Name, Month, DNA_or_cDNA, Fox_Aspect, Snowpack_Layer))
ps.final.measures = estimate_richness(ps.final, measures = c("Observed","Shannon","Simpson")) %>%
rownames_to_column(var = "Sample_ID") %>%
mutate(Sample_ID = gsub("\\.", "-", Sample_ID))
full_join(ps.final.measures, mini_metadata,
by="Sample_ID") %>%
kable(digits = 4) %>%
kable_styling(full_width = T, bootstrap_options = c("striped", "hover")) %>%
row_spec(0, bold=TRUE)| Sample_ID | Observed | Shannon | Simpson | Sample_Name | Month | DNA_or_cDNA | Fox_Aspect | Snowpack_Layer |
|---|---|---|---|---|---|---|---|---|
| A03-cDNA-1 | 286 | 0.7114 | 0.2103 | April_N_NE_Basal_R | April | cDNA | N_NE | BASAL |
| A03-DNA-1 | 138 | 3.2774 | 0.9407 | June_N_NE_Top20_D | June | DNA | N_NE | TOP |
| A04-cDNA-2 | 163 | 2.0361 | 0.7852 | April_NW_AWS_Top_R | April | cDNA | NW_AWS | TOP |
| A04-DNA-2 | 75 | 1.6409 | 0.6083 | April_ NW_AWS_ Basal_D | April | DNA | NW_AWS | BASAL |
| A05-cDNA-3 | 515 | 2.2067 | 0.7048 | June_NW_AWS_Basal_R | June | cDNA | NW_AWS | BASAL |
| A05-DNA-3 | 1390 | 5.1230 | 0.9813 | 30 July_ N_NE_ Gl ice_D | Late July | DNA | N_NE | GL ICE |
| A06-cDNA-4 | 125 | 1.3510 | 0.6058 | June_N_NE_Basal_R | June | cDNA | N_NE | BASAL |
| A06-DNA-4 | 210 | 3.9337 | 0.9650 | June_ SWS1SE_ Top20_D | June | DNA | SWS1SE | TOP |
| A07-cDNA-5 | 1635 | 4.6591 | 0.9488 | 9 July_NW_AWS_Mid_R | Early July | cDNA | NW_AWS | MID |
| A07-DNA-5 | 609 | 3.3602 | 0.8726 | 9 July_ NW_AWS_ Sup ice_D | Early July | DNA | NW_AWS | SUP ICE |
| A08-DNA-6 | 563 | 4.0141 | 0.8786 | April_ NW_AWS_ Mid_D | April | DNA | NW_AWS | MID |
| A09-cDNA-7 | 1875 | 5.5902 | 0.9807 | 30 July_NW_AWS_Top20_R | Late July | cDNA | NW_AWS | TOP |
| A10-cDNA-8 | 842 | 4.7339 | 0.9520 | 30 July_N_NE_Gl ice_R | Late July | cDNA | N_NE | GL ICE |
| A10-DNA-8 | 179 | 3.3916 | 0.9493 | April_NW_AWS_Top20_D | April | DNA | NW_AWS | TOP |
| A11-cDNA-9 | 3 | 0.2518 | 0.1121 | April_NW_AWS_Mid_R | April | cDNA | NW_AWS | MID |
| A11-DNA-9 | 1819 | 5.3121 | 0.9764 | 9 July_SWS1SE_Gl ice 1_D | Early July | DNA | SWS1SE | GL ICE |
| A12-cDNA-10 | 54 | 1.4929 | 0.6133 | June_NW_AWS_Top20_R | June | cDNA | NW_AWS | TOP |
| A12-DNA-10 | 2719 | 5.7179 | 0.9837 | 30 July_N_NE_Top20_D | Late July | DNA | N_NE | TOP |
| B01-cDNA-11 | 453 | 4.0612 | 0.9622 | 30 July_NW_AWS_Sup ice_R | Late July | cDNA | NW_AWS | SUP ICE |
| B01-DNA-11 | 187 | 3.5277 | 0.9546 | June_NW_AWS_Mid_D | June | DNA | NW_AWS | MID |
| B02-cDNA-12 | 733 | 4.5851 | 0.9683 | 30 July_N_NE_Top20_R | Late July | cDNA | N_NE | TOP |
| B02-DNA-12 | 310 | 2.1151 | 0.7154 | 9 July_N_NE_Mid_D | Early July | DNA | N_NE | MID |
| B03-cDNA-13 | 38 | 1.1241 | 0.4091 | June_SWS1SE_Basal_R | June | cDNA | SWS1SE | BASAL |
| B03-DNA-13 | 457 | 4.2859 | 0.9559 | April_SWS1SE_Mid_D | April | DNA | SWS1SE | MID |
| B04-cDNA-14 | 67 | 1.1867 | 0.4139 | 9 July_NW_AWS_Top20_R | Early July | cDNA | NW_AWS | TOP |
| B05-cDNA-15 | 1307 | 5.7881 | 0.9881 | 9 July_NW_AWS_Gl ice 1_R | Early July | cDNA | NW_AWS | GL ICE |
| B06-cDNA-16 | 134 | 1.4982 | 0.6595 | April_SWS1SE_Top20_R | April | cDNA | SWS1SE | TOP |
| B06-DNA-2 | 156 | 0.6586 | 0.2271 | April_SWS1SE_Basal_D | April | DNA | SWS1SE | BASAL |
| B07-DNA-3 | 1302 | 5.3971 | 0.9910 | June_NW AWS_Top20_D | June | DNA | NW_AWS | TOP |
| B08-cDNA-18 | 21 | 0.9250 | 0.5154 | June_NW_AWS_Mid_R | June | cDNA | NW_AWS | MID |
| B08-DNA-4 | 986 | 5.2137 | 0.9805 | 30 July_SWS1SE_Gl ice_D | Late July | DNA | SWS1SE | GL ICE |
| B09-DNA-5 | 105 | 2.9238 | 0.9150 | June_N_NE_Mid_D | June | DNA | N_NE | MID |
| B10-cDNA-21 | 400 | 4.5700 | 0.9771 | 9 July_SWS1SE_Sup ice_R | Early July | cDNA | SWS1SE | SUP ICE |
| B10-DNA-6 | 377 | 4.1547 | 0.9598 | June_SWS1SE_Basal_D | June | DNA | SWS1SE | BASAL |
| B12-cDNA-23 | 29 | 0.9625 | 0.3511 | 30 July_SWS1SE_Slush_R | Late July | cDNA | SWS1SE | TOP |
| B12-DNA-8 | 1998 | 5.3382 | 0.9817 | 9 July_ NW_AWS_Gl ice 1_D | Early July | DNA | NW_AWS | GL ICE |
| C01-DNA-9 | 165 | 3.3370 | 0.9306 | June_SWS1SE_Mid_D | June | DNA | SWS1SE | MID |
| C02-cDNA-26 | 907 | 4.1753 | 0.9419 | 30 July_SWS1SE_Gl ice_R | Late July | cDNA | SWS1SE | GL ICE |
| C03-cDNA-27 | 537 | 3.3393 | 0.8840 | 9 July_N_NE_Sup ice_R | Early July | cDNA | N_NE | SUP ICE |
| C03-DNA-B2-1 | 677 | 5.1429 | 0.9887 | June_NW_AWS_Basal_D | June | DNA | NW_AWS | BASAL |
| C04-cDNA-28 | 307 | 3.3747 | 0.8749 | 30 July_NW_AWS_Gl ice_R | Late July | cDNA | NW_AWS | GL ICE |
| C04-DNA-2 | 715 | 4.2502 | 0.9430 | 30 July_NW_AWS_Gl ice_D | Late July | DNA | NW_AWS | GL ICE |
| C05-DNA-3 | 1739 | 5.2347 | 0.9838 | 30 July_NW_AWS_Sup ice_D | Late July | DNA | NW_AWS | SUP ICE |
| C06-DNA-4 | 236 | 3.4402 | 0.9006 | April_SWS1SE_Top20_D | April | DNA | SWS1SE | TOP |
| C07-DNA-5 | 2789 | 5.9749 | 0.9894 | 9 July_NW_AWS_Mid_D | Early July | DNA | NW_AWS | MID |
| C08-DNA-6 | 307 | 4.5001 | 0.9837 | June_N_NE_Basal_D | June | DNA | N_NE | BASAL |
| C09-cDNA-33 | 36 | 1.3284 | 0.5802 | June_N_NE_Top20_R | June | cDNA | N_NE | TOP |
| C09-DNA-7 | 148 | 3.8243 | 0.9690 | 9 July_N_NE_Gl ice_D | Early July | DNA | N_NE | GL ICE |
| C10-DNA-8 | 383 | 3.2676 | 0.8801 | 9 July_ SWS1SE_Mid_D | Early July | DNA | SWS1SE | MID |
| C11-cDNA-35 | 65 | 1.5747 | 0.6657 | 9 July_SWS1SE_Mid_R | Early July | cDNA | SWS1SE | MID |
| C12-cDNA-36 | 19 | 0.4213 | 0.1462 | 9 July_N_NE_Gl ice_R | Early July | cDNA | N_NE | GL ICE |
| C12-DNA-10 | 396 | 2.2692 | 0.7492 | 9 July_N_NE_Sup ice_D | Early July | DNA | N_NE | SUP ICE |
| D01-DNA-B4-1 | 147 | 0.6873 | 0.2422 | 30 July_SWS1SE_Slush_D | Late July | DNA | SWS1SE | TOP |
| D02-cDNA-38 | 128 | 2.1030 | 0.7780 | April_SWS1SE_Mid_R | April | cDNA | SWS1SE | MID |
| D03-cDNA-39 | 266 | 3.4007 | 0.9035 | April_N_NE_Mid_R | April | cDNA | N_NE | MID |
| D04-cDNA-40 | 79 | 1.5083 | 0.6114 | June_SWS1SE_Top20_R | June | cDNA | SWS1SE | TOP |
| D05-cDNA-41 | 154 | 1.8281 | 0.7377 | June_N_NE_Mid_R | June | cDNA | N_NE | MID |
| D05-DNA-5 | 90 | 0.6634 | 0.2169 | 9 July_N_NE_Top20_D | Early July | DNA | N_NE | TOP |
| D06-cDNA-42 | 841 | 3.6042 | 0.9024 | 9 July_NW_AWS_Sup ice_R | Early July | cDNA | NW_AWS | SUP ICE |
| D06-DNA-6 | 2292 | 5.5649 | 0.9810 | 9 July_SWS1SE_Sup ice_D | Early July | DNA | SWS1SE | SUP ICE |
| D07-cDNA-43 | 1444 | 4.5026 | 0.9489 | 9 July_SWS1SE_Gl ice 2_R | Early July | cDNA | SWS1SE | GL ICE |
| D07-DNA-7 | 467 | 1.8746 | 0.5308 | 9 July_NW_AWS_Top20_D | Early July | DNA | NW_AWS | TOP |
| D08-cDNA-44 | 120 | 2.6011 | 0.8428 | 9 July_N_NE_Mid_R | Early July | cDNA | N_NE | MID |
| D08-DNA-8 | 1045 | 4.1510 | 0.9245 | 30 July_NW_AWS_Top20_D | Late July | DNA | NW_AWS | TOP |
| D09-DNA-9 | 2986 | 5.9104 | 0.9895 | 9 July_SWS1SE_Top20_D | Early July | DNA | SWS1SE | TOP |
| B05-DNA-B1-1 | NA | NA | NA | April_N_NE_Top20_D | April | DNA | N_NE_D | TOP |
| C02-DNA-10 | NA | NA | NA | April_N_NE_Mid_D | April | DNA | N_NE_D | MID |
| C11-DNA-9 | NA | NA | NA | April_N_NE_Basal_D | April | DNA | N_NE_D | BASAL |
| D02-DNA-2 | NA | NA | NA | 9 July_NW_AWS_Gl ice 2_D | Early July | DNA | NW_AWS | GL ICE |
| B11-DNA-7 | NA | NA | NA | 9 July_SWS1SE_Gl ice 2_D | Early July | DNA | SWS1SE | GL ICE |
| A01-DNA-MM | NA | NA | NA | MM_D | PCR_Control | DNA | Control | Control |
| A02-DNA-MM-H20 | NA | NA | NA | MM+H2O_D | PCR_Control | DNA | Control | Control |
| A09-DNA-7 | NA | NA | NA | St MQ_D | Control | DNA | Control | Control |
| B04-DNA-14 | NA | NA | NA | St MQ Unis_D | Control | DNA | Control | Control |
| D03-DNA-3 | NA | NA | NA | SLS_D | Control | DNA | Control | Control |
| D04-DNA-4 | NA | NA | NA | ST Control_D | Control | DNA | Control | Control |
| D10-DNA-B5-1 | NA | NA | NA | St MQ_D | Control | DNA | Control | Control |
| D11-DNA-2 | NA | NA | NA | SLS St_D | Control | DNA | Control | Control |
| D12-DNA-3 | NA | NA | NA | St Only_D | Control | DNA | Control | Control |
| E01-DNA-4 | NA | NA | NA | St MQ2_D | Control | DNA | Control | Control |
| E02-DNA-MM | NA | NA | NA | MM_D | PCR_Control | DNA | Control | Control |
| E03-DNA-MM-H20 | NA | NA | NA | MM+H20_D | PCR_Control | DNA | Control | Control |
| C07-cDNA-31 | NA | NA | NA | April_NW_AWS_Basal_R | April | cDNA | NW_AWS | BASAL |
| B07-cDNA-17 | NA | NA | NA | April_N_NE_Top20_R | April | cDNA | N_NE | TOP |
| C08-cDNA-32 | NA | NA | NA | April_SWS1SE_Basal_R | April | cDNA | SWS1SE | BASAL |
| B09-cDNA-19 | NA | NA | NA | June_SWS1SE_Mid_R | June | cDNA | SWS1SE | MID |
| C10-cDNA-34 | NA | NA | NA | 9 July_NW_AWS_Gl ice 2_R | Early July | cDNA | NW_AWS | GL ICE |
| B11-cDNA-22 | NA | NA | NA | 9 July_ N_NE_Top20_R | Early July | cDNA | N_NE | TOP |
| A08-cDNA-6 | NA | NA | NA | 9 July_SWS1SE_Gl ice_R | Early July | cDNA | SWS1SE | GL ICE |
| A01-cDNA-MM | NA | NA | NA | MM_R | PCR_Control | cDNA | Control | Control |
| A02-cDNA-MM-H20 | NA | NA | NA | MM+H20_R | PCR_Control | cDNA | Control | Control |
| C01-cDNA-25 | NA | NA | NA | St MQ Unis_R | Control | cDNA | Control | Control |
| C05-cDNA-29 | NA | NA | NA | SLS St_R | Control | cDNA | Control | Control |
| C06-cDNA-30 | NA | NA | NA | St MQ2_R | Control | cDNA | Control | Control |
| D01-cDNA-37 | NA | NA | NA | St Only_R | Control | cDNA | Control | Control |
| D09-cDNA-45 | NA | NA | NA | SLS St_R | Control | cDNA | Control | Control |
| D10-cDNA-46 | NA | NA | NA | St MQ_R | Control | cDNA | Control | Control |
| D11-cDNA-47 | NA | NA | NA | St Only_R | Control | cDNA | Control | Control |
| D12-cDNA-MM | NA | NA | NA | MM_R | PCR_Control | cDNA | Control | Control |
| E01-cDNA-MM-H20 | NA | NA | NA | MM+H20_R | PCR_Control | cDNA | Control | Control |
plot_richness(ps.final,
x="Month", color="Month",
measures=c("Observed","Shannon","Simpson")) +
geom_violin(alpha = .5) +
geom_jitter(shape=16, position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = month_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
labs(color="Month") +
theme_bw() +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))plot_richness(subset_samples(ps.final,
Snowpack_Layer != "GL ICE" &
Sample_Name != "30 July_SWS1SE_Slush_R" &
Sample_Name != "30 July_SWS1SE_Slush_D" &
!(Sample_Name %in% gl_ice_spls)),
x="Month", color="Month",
measures=c("Observed","Shannon","Simpson")) +
geom_violin(alpha = .5) +
geom_jitter(shape=16, position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = month_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
labs(color="Month") +
theme_bw() +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))mini_month_comparisons = list(c("April", "June"),
c("April", "June"))
plot_richness(subset_samples(ps.final, Month != "Early July" & Month != "Late July"),
x="Month", color="Month",
measures=c("Observed","Shannon","Simpson")) +
geom_violin(alpha = .5) +
geom_jitter(shape=16, position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = mini_month_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
labs(color="Month") +
theme_bw() +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))#now without stats, only shapes
plot_richness(ps.final,
x="Fox_Aspect", color="Fox_Aspect",
measures=c("Observed","Shannon","Simpson")) +
geom_violin() +
geom_jitter(aes(shape=Month), size=5, alpha=0.7,
position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
labs(color="Fox Aspect") +
scale_shape_manual(values=nmds_shapes) +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))plot_richness(ps.final,
x="Fox_Aspect", color="Fox_Aspect",
measures=c("Observed","Shannon","Simpson")) +
geom_violin(alpha=0.5) +
geom_jitter(shape=16, position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = site_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
labs(color="Fox Aspect") +
theme_bw() +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))plot_richness(subset_samples(ps.final,
Snowpack_Layer != "GL ICE" &
Sample_Name != "30 July_SWS1SE_Slush_R" &
Sample_Name != "30 July_SWS1SE_Slush_D" &
!(Sample_Name %in% gl_ice_spls)),
x="Fox_Aspect", color="Fox_Aspect",
measures=c("Observed","Shannon","Simpson")) +
geom_violin(alpha=0.5) +
geom_jitter(shape=16, position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = site_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
labs(color="Fox Aspect") +
theme_bw() +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))plot_richness(subset_samples(ps.final, Snowpack_Layer != "GL ICE" &
Month != "Early July" & Month != "Late July"),
x="Fox_Aspect", color="Fox_Aspect",
measures=c("Observed","Shannon","Simpson")) +
geom_violin(alpha = .5) +
geom_jitter(shape=16, position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = mini_month_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
labs(color="Fox_Aspect") +
theme_bw() +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
plot_richness(subset_samples(ps.final, Month == c("Early July", "Late July")),
x="Fox_Aspect", color="Fox_Aspect",
measures=c("Observed","Shannon","Simpson")) +
geom_violin(alpha = .5) +
geom_jitter(shape=16, position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = mini_month_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
labs(color="Fox_Aspect") +
theme_bw() +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))## Warning in `==.default`(Month, c("Early July", "Late July")): longer object
## length is not a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
plot_richness(subset_samples(ps.final, Month == c("Early July", "Late July"),
Snowpack_Layer != "GL ICE" &
Sample_Name != "30 July_SWS1SE_Slush_R" &
Sample_Name != "30 July_SWS1SE_Slush_D" &
!(Sample_Name %in% gl_ice_spls)),
x="Fox_Aspect", color="Fox_Aspect",
measures=c("Observed","Shannon","Simpson")) +
geom_violin(alpha = .5) +
geom_jitter(shape=16, position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = mini_month_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
labs(color="Fox_Aspect") +
theme_bw() +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))## Warning in `==.default`(Month, c("Early July", "Late July")): longer object
## length is not a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
plot_richness(subset_samples(ps.final, Snowpack_Layer != "GL ICE" &
Month == "Early July" & !(Sample_Name %in% gl_ice_spls)),
x="Fox_Aspect", color="Fox_Aspect",
measures=c("Observed","Shannon","Simpson")) +
geom_violin(alpha = .5) +
geom_jitter(shape=16, position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = mini_month_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
labs(color="Fox_Aspect") +
theme_bw() +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
plot_richness(subset_samples(ps.final, Snowpack_Layer != "GL ICE" & Month == "Late July" &
Sample_Name != "30 July_SWS1SE_Slush_R" &
Sample_Name != "30 July_SWS1SE_Slush_D" &
!(Sample_Name %in% gl_ice_spls)),
x="Fox_Aspect", color="Fox_Aspect",
measures=c("Observed","Shannon","Simpson")) +
geom_violin(alpha = .5) +
geom_jitter(shape=16, position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = mini_month_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
labs(color="Fox_Aspect") +
theme_bw() +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
#now without stats, only shapes
plot_richness(ps.final,
x="Snowpack_Layer", color="Snowpack_Layer",
measures=c("Observed","Shannon","Simpson")) +
geom_violin() +
geom_jitter(aes(shape=Month), size=5, alpha=0.7,
position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
labs(color="Snowpack Layer") +
scale_shape_manual(values=nmds_shapes) +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))plot_richness(ps.final,
x="Snowpack_Layer", color="Snowpack_Layer",
measures=c("Observed","Shannon","Simpson")) +
geom_violin(alpha = .5) +
geom_jitter(position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = slayer_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
labs(color="Snowpack Layer") +
theme_bw() +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))plot_richness(subset_samples(ps.final, Month = c("Early July", "Late July")),
x="Snowpack_Layer", color="Snowpack_Layer",
measures=c("Observed","Shannon","Simpson")) +
geom_violin(alpha = .5) +
geom_jitter(position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = slayer_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
labs(color="Snowpack Layer") +
theme_bw() +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))# ggsave("july_alpha.svg", width = 12, height = 8)
# ggsave("july_alpha.tiff", width = 12, height = 8)plot_richness(subset_samples(ps.final, Month == c("Early July", "Late July") &
Snowpack_Layer != "GL ICE" &
!(Sample_Name %in% gl_ice_spls) &
Sample_Name != c("30 July_SWS1SE_Slush_R",
"30 July_SWS1SE_Slush_D")),
x="Snowpack_Layer", color="Snowpack_Layer",
measures=c("Observed","Shannon","Simpson")) +
geom_violin(alpha = .5) +
geom_jitter(position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = slayer_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
labs(color="Snowpack Layer") +
theme_bw() +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))## Warning in `==.default`(Month, c("Early July", "Late July")): longer object
## length is not a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length
## Warning in `!=.default`(Sample_Name, c("30 July_SWS1SE_Slush_R", "30
## July_SWS1SE_Slush_D")): longer object length is not a multiple of shorter object
## length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length
## Warning: Computation failed in `stat_signif()`:
## not enough 'x' observations
## Warning: Computation failed in `stat_signif()`:
## not enough 'x' observations
## Warning: Computation failed in `stat_signif()`:
## not enough 'x' observations
## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations
## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations
## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations
# ggsave("july_alpha.svg", width = 12, height = 8)
# ggsave("july_alpha.tiff", width = 12, height = 8)plot_richness(subset_samples(ps.final, Snowpack_Layer != "GL ICE" &
Month == "Early July" & !(Sample_Name %in% gl_ice_spls)),
x="Snowpack_Layer", color="Snowpack_Layer",
measures=c("Observed","Shannon","Simpson")) +
geom_violin(alpha = .5) +
geom_jitter(position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = slayer_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
labs(color="Snowpack Layer") +
theme_bw() +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## missing value where TRUE/FALSE needed
## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations
## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations
## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations
# ggsave("july_alpha.svg", width = 12, height = 8)
# ggsave("july_alpha.tiff", width = 12, height = 8)plot_richness(subset_samples(ps.final, Snowpack_Layer != "GL ICE" &
Sample_Name != "30 July_SWS1SE_Slush_R" &
Sample_Name != "30 July_SWS1SE_Slush_D" &
Month == "Late July" & !(Sample_Name %in% gl_ice_spls)),
x="Snowpack_Layer", color="Snowpack_Layer",
measures=c("Observed","Shannon","Simpson")) +
geom_violin(alpha = .5) +
geom_jitter(position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = slayer_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
labs(color="Snowpack Layer") +
theme_bw() +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12, face = "bold"),
axis.text.y = element_text(face = "bold"),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`:
## replacement has 1 row, data has 0
## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations
## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations
## Warning: Computation failed in `stat_signif()`:
## not enough 'y' observations
## Warning: Computation failed in `stat_signif()`:
## not enough 'x' observations
## Warning: Computation failed in `stat_signif()`:
## not enough 'x' observations
## Warning: Computation failed in `stat_signif()`:
## not enough 'x' observations
# ggsave("july_alpha.svg", width = 12, height = 8)
# ggsave("july_alpha.tiff", width = 12, height = 8)ps.final.rare = rarefy_even_depth(ps.final,
sample.size = min(sample_sums(ps.final)),
rngseed = TRUE, replace = TRUE,
trimOTUs = TRUE, verbose = TRUE)## `set.seed(TRUE)` was used to initialize repeatable random subsampling.
## Please record this for your records so others can reproduce.
## Try `set.seed(TRUE); .Random.seed` for the full vector
## ...
## 25951OTUs were removed because they are no longer
## present in any sample after random subsampling
## ...
plot_richness(ps.final.rare,
x="Sample_ID", color="Sample_ID",
measures=c("Observed","Shannon","Simpson")) +
geom_jitter(shape=16, position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
labs(color="Sample_ID") +
guides(color=guide_legend(ncol=2)) +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))plot_richness(ps.final.rare,
x="Month", color="Month",
measures=c("Observed","Shannon","Simpson")) +
geom_violin() +
geom_jitter(shape=16, position=position_jitter(0.2)) +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = month_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
labs(color="Month") +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))plot_richness(ps.final.rare,
x="Fox_Aspect", color="Fox_Aspect",
measures=c("Observed","Shannon","Simpson")) +
geom_violin() +
geom_jitter(shape=16, position=position_jitter(0.2)) +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = site_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
labs(color="Fox Aspect") +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))#now without stats, only shapes
plot_richness(ps.final.rare,
x="Fox_Aspect", color="Fox_Aspect",
measures=c("Observed","Shannon","Simpson")) +
geom_violin() +
geom_jitter(aes(shape=Month), size=5, alpha=0.7,
position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
labs(color="Fox Aspect") +
scale_shape_manual(values=nmds_shapes) +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))plot_richness(ps.final.rare,
x="Snowpack_Layer", color="Snowpack_Layer",
measures=c("Observed","Shannon","Simpson")) +
geom_violin() +
geom_jitter(position=position_jitter(0.2)) +
# stat_compare_means(comparisons = my_comparisons)+
stat_compare_means(comparisons = slayer_comparisons,
method = "t.test",
var.equal="welch",
label.x.npc = "right",
label.y.npc = "top",
symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05,0.01, 1),
symbols = c("p<0.0001", "p<0.001", "p<0.01", "p<0.05","p<0.1", "ns"))) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
labs(color="Snowpack Layer") +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))#now without stats, only shapes
plot_richness(ps.final.rare,
x="Snowpack_Layer", color="Snowpack_Layer",
measures=c("Observed","Shannon","Simpson")) +
geom_violin() +
geom_jitter(aes(shape=Month), size=5, alpha=0.7,
position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
labs(color="Snowpack Layer") +
scale_shape_manual(values=nmds_shapes) +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))plot_richness(ps.final.rare,
x="Snowpack_Layer", color="Snowpack_Layer",
measures=c("Observed","Shannon","Simpson")) +
geom_violin() +
geom_jitter(aes(shape=Fox_Aspect), size=5, alpha=0.7,
position=position_jitter(0.2)) +
facet_wrap(DNA_or_cDNA~variable, scales = "free") +
labs(color="Snowpack Layer") +
scale_shape_manual(values=nmds_shapes) +
theme(axis.text.x = element_text(angle = 20,
vjust = 1,
hjust = 1,
size = 12),
axis.title = element_blank(),
strip.text = element_text(size = 12, margin = margin()))